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Abstract 

The stability of the axisymmetric solitary waves of the Gross-Pitaevskii (GP) equa- 
tion is investigated. The Implicitly Restarted Arnoldi Method for banded matrices with 
shift-invert was used to solve the linearised spectral stability problem. The rarefaction 
solitary waves on the upper branch of the Jones-Roberts dispersion curve are shown to 
be unstable to axisymmetric infinitesimal perturbations, whereas the solitary waves on 
the lower branch and all two-dimensional solitary waves are linearly stable. The growth 
rates of the instabilities on the upper branch are so small that an arbitrarily specified 
initial perturbation of a rarefaction wave at first usually evolves towards the upper 
branch as it acoustically radiates away its excess energy. This is demonstrated through 
numerical integrations of the GP equation starting from an initial state consisting of an 
unstable rarefaction wave and random non-axisymmetric noise. The resulting solution 
evolves towards, and remains for a significant time in the vicinity of, an unperturbed 
unstable rarefaction wave. It is shown however that, ultimately (or for an initial state 
extremely close to the upper branch), the solution evolves onto the lower branch or is 
completely dissipated as sound. 

Pacs: 03.75.Lm, 05.45.-a, 67.40.Vs, 67.57.De 

1 Introduction 

Theoretical investigation of the structure, energy, dynamics, and stability of vortices in 
Bose-Einstein condensates (BEC) has received increased attention since Bose-Einstein con- 
densation was achieved in trapped alkali-metal gases at ultralow temperatures; for a com- 
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prehensive review, see for example [Q. Quantised vorticity is an intriguing feature of su- 
perfluidity, and much effort was therefore devoted to manipulating vortices and observing 
their dynamics. This provided a valuable tool for elucidating the physics of many-particle 
systems and relating it to the quantitative predictions of thermal field theories. The con- 
densates of alkali vapours are pure and dilute, so that the Gross-Pitaevskii (GP) model, 
which is the so-called 'mean-field limit' of quantum field theories, gives a precise description 
of their dynamics at low temperatures. The structure and properties of quantised vortices 
were therefore studied experimentally and then refined, both analytically and numerically, 
by using the GP equation. Conversely, the predictions on vortex properties based on the 
GP model were later confirmed experimentally. Besides vortices, only one other localised 
disturbance has been detected experimentally: the one-dimensional solitary wave (the 'dark 
soliton') was created from a density depletion (phase imprinting) in a BEC of sodium atoms 

0- 

The focus of this report is on a different class of solitary waves that should exist in 
a condensate: these are the vorticity-free axisymmetric disturbances whose existence was 
predicted by Jones and Roberts [3] on the basis of numerical integrations of 

-2id t v = vV + (i-|Vf)V>, (i) 

where dt = d/dt. In this dimensionless form of the GP equation, the unit of length 
is the healing length, the speed of sound is c = l/\/2, and the density at infinity is 
Poo = iV'ool 2 = 1- In contrast with the vortices and vortex rings that were extensively 
studied in superfluid helium systems long before the experimental realisation of Bose- 
Einstein condensation, the rarefaction solitary waves do not exist in superfluid helium 
as they move faster than the Landau critical velocity. Therefore, the existence of the 
rarefaction solitary waves is a unique characteristic of BEC. 

Jones and Roberts found all solitary wave solutions of the GP equation in two dimen- 
sions (2D) and three dimensions (3D). In a momentum energy (pE) plot, the 3D sequence, 
which we call the 'JR dispersion curve', has two branches meeting at a cusp where p and 
E simultaneously assume their minimum values, p c and E c . Here 

p = jJ\(r-i)d z ^-^-i)d z r]dv, (2) 

E = l -J\V^dV+ I|(l-|Vf) 2 ^, (3) 

where z is the direction in which the wave propagates and the direction about which it 
is axisymmetric. For each p in excess of the minimum p c , two values of E are possible, 
and £^ooasp->ooon each branch. On the lower (energy) branch the solutions are 
asymptotic to large circular vortex rings. As p and E decrease from infinity on this branch, 
the solutions begin to lose their similarity to vortex rings. Eventually, for a momentum po 
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Figure 1: (Color online) The JR dispersion curve for the family of the axisymmetric solitary 
wave solutions. The part of the curve that corresponds to vortex rings is shown in grey 
(red). Three rarefaction waves that are considered in the text are indicated by circles. 
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slightly in excess of p c , they lose their vorticity (tp loses its zero), and thereafter the solutions 
may better be described as 'rarefaction waves'. The upper branch solutions consist entirely 
of these waves and, as p — > oo, they asymptotically approach the rational soliton solution 
of the Kadomtsev-Petviashvili Type I (KPI) equation. Fig. ^ shows the two branches, the 
cusp, and the three rarefaction waves that will provide the principal examples below. 

The JR sequence can be uniquely characterised in other ways, for example by the 
velocity, U, of the wave or by min(n), the minimum of the real part of ip = u + iv. Fig. [2] 
shows min(u) as a function of U for the entire family of JR solutions. In the limit U — > of 
large vortex rings, min(u) — ► —1; in the opposite limit, U — > c, of large rarefaction waves, 
min(u) — > +1. Between these extremes, the case min(it) = deserves special mention, 
as it separates the rarefaction waves, which do not possess vorticity, from the vortex- type 
solutions which do. In this case, ip vanishes at a single point, so that this solution might 
appropriately be termed a 'point defect'; its velocity is approximately U ~ 0.62. The cusp 
(U ~ 0.65) in the pE— plot arises because U = dE/dp = E'(U)/p'(U), so that extrema of 
p are simultaneously extrema of E. 

It was suggested in :.3. as well as in its more detailed sequel jl] that every solitary 
wave on the upper branch is unstable, since it is energetically favourable for it to 'collapse' 
onto the lower branch of smaller energy at the same momentum. A Derrick-type argument 
[3] was used in which neighbouring axisymmetric states having the same p as the upper 
branch solution to (1) were shown to have a smaller E. It would follow that the solution 
is unstable provided that p and E are the only quantities conserved by (1) in 3D, but this 
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Figure 2: The minimum of the real part of the wavefunction, Re{ip) = u, for the JR solitary 
waves as a function of velocity U. The open circle marks the point defect; closed circles 
show computed points through which the continuous line is drawn. 




has never been demonstrated. Moreover, rarefaction waves on the upper branch of the JR 
dispersion curve are seen in numerical simulations: they evolve from a density depletion of 
a condensate JU] and appear during condensate formation from a strongly non-equilibrated 
Bose gas |14j . They may exist for a considerable time before they either collapse onto the 
lower branch of the JR dispersion curve or lose their energy to sound waves and disappear. 

Evidently there is a paradox. On the one hand, if the assumption that p and E are the 
only conserved quantities in 3D is correct and the upper branch solutions are intrinsically 
unstable, how could they even form? On the other hand, if the assumption is false and the 
upper branch solutions are stable, why do they ultimately disappear? Is it because they 
are metastable or because they lose energy and momentum through collisions with other 
waves in the system? 

The goal of this paper is to resolve this paradox. In Section 2, we demonstrate that 
the upper branch solitary waves are unstable. In Section 3.1, we show that, nevertheless, a 
perturbed solitary wave on the upper branch will generally evolve towards an unperturbed 
state on that branch and remain in its vicinity for a long period. Eventually, however, as 
illustrated in Section 3.2, it collapses onto the lower branch or disappears entirely. Another 
facet of such evolutionary processes is studied in Section 3.3, where it is shown how two 
rarefaction waves from the upper branch that move in the same direction can create a 
vortex ring. In Section 3.4 we demonstrate how unstable rarefaction waves can form from 
the evolution of a density depletion. We conclude in Section 4 with a brief summary of our 
findings. 



4 



2 Linear stability of the rarefaction solitary waves 

2.1 The eigenvalue problem 

In a reference frame moving with velocity U in the z— direction, the GP equation Q is 

-2idt1> + 2iUd z ip = V 2 V + (i - iVf (4) 

solutions to which must obey 

tp — > 1, for r = |x| — > oo . (5) 

The "basic solution" to @ and ((SJ) is the solitary wave, ipo (= no + Wq), for which 

2iZ7S,^o = vVo + ll-lV'ol 2 )^, (6) 
■00 — ^ 1; f° r r — ► oo . (7) 

In this Section, we seek to determine the fate of small perturbations to tpQ. We write 

^(x,t) = Vo(x) + ^(x,t). (8) 

On substituting Q into (jlj) and linearising with respect to ip, we find that 

- 2idS + 2iUd z $ = V 2 ^ + (1 - 2|Vo| 2 )^ - ^* , (9) 
— ► 0, for r — > oo , (10) 

where * stands for complex conjugation. Further discussion of (jTU)) is postponed to §2.4. 
The objective now is the general solution of the initial value problem, i.e., the solution of Q 
and (|l(Jj) . starting from an arbitrarily specified initial state. Following a well-trodden path 
in stability theory, we suppose that the solution can be expressed as a linear combination 
of a complete set of normal mode solutions. Alternative approaches to the initial value 
problem will be described in §3. 

2.2 Real growth rates 

Equation © admits normal mode solutions in which 

^(x,t) =Pi(x)e CTi , ^(x,t) =^(x)e CTt , (11) 

and 

Lpi ~ 4>IpI = -2i<rpi , L*p* 2 - % 2 p x = 2\ap* 2 , (12) 

where 

L = V 2 + 1 - 2|V> | 2 - 2iUd z . (13) 
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If the growth rate a is real, we may take 

Pi = P2 = P = u + iv (say) , 



(14) 



where u and v are real, both being proportional to cosm</> or sinm0 where m is an in- 
teger and (s, 4>, z) are cylindrical coordinates. We shall be primarily interested in the 
axisymmetric case m = 0. Equations (|12|) require 



Lp - t/>o£»* 



-2icrp, 



or equivalently 



V 2 u + 2Ud z v + (1 - 3uq - v%)u - 2u v v 
V 2 v — 2Ud z u + (1 — Uq — Svfyv — 2uqVqu 
u,v — ► 0, for r 



2av , 
-2au , 



oo . 



(15) 



(16) 
(17) 
(18) 



where ipo = uo + • 

The numerical solution of (|16|) - l)18[) . was achieved by first making the change of 
independent variables suggested by [2|: z' = z and s' = s-y/(l — 2[/ 2 ). We mapped the 
infinite domain onto the box (0, ^vr) x (— hn, ^vr) using the transformation 



z = tan l [Cz), 



(19) 



with a similar transformation for s'; here C (~ 0.4 — 0.5) is a constant chosen at our 
convenience. We also write ip = + 1, so that the boundary conditions at infinity are 
4» = 0. 

We used the Newton- Raphson iteration procedure (using a banded matrix linear solver 
based on the bi-conjugate gradient stabilised iterative method with preconditioning) in the 
frame of reference moving with velocity U to find solitary waves from solutions of 



2iURcos 2 zd^ 
and 



- (*o + % + l^o| 2 )(l + *o) 
0, as s, z — ► oo, 



(20) 



where 



C 



R 2 (l 



2U 2 )cos s s 



cos s 



ds 2 



2 sins 



sm s 



d_ 

ds 



+ R 2 cos 3 £ 



cos z 



d^ 2 



^ d 
2 sm z 

oz 



(21) 



We also used (j!9|) to transform (|16|) - (|18[). The resulting matrix equation was then 
discretised and gave the classic eigenvalue problem, Ax. = erx, for a. Eigenvalues can often 
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be found successively by the Arnoldi/Lanczos process, in which the extreme eigenvalues are 
found first and the remaining eigenvalues are found successively by the 'shift and invert' 
technique, in which the original eigenvalue is shifted by A (say), and A — XI is inverted to 
transform the equation to (A — A/) _1 x = fix, where [i = \j{p — A), so that the original 
eigenvalue is easily recovered as + A. We employed the ARPACK collection of 
subroutines, that uses the Implicitly Shifted QR technique that is suitable for large scale 
problems. The algorithm is designed to compute a few (k) eigenvalues with user specified 
features such as those of largest real part or largest magnitude. Storage requirements 
are on the order of n * k locations, where n is a number of columns of the matrix A. A 
set of Schur basis vectors for the desired £;-dimensional eigenspace is computed which is 
numerically orthogonal to working precision. 

The basic solution was found with a maximum resolution of 250x200, making the real 
matrix generated by (|16|) - (|17j) of order 10 5 x 10 5 . The resolution of the basic solution 
was halved to check the accuracy of the eigenvalues found. 

We found that a 2 is real for all solutions of (|16|) - (|18|) , To explain this, we first observe 
that we are interested only in perturbations that preserve the mass of the state ip$ in the 
sense that, if use (JHJ) to evaluate the mass flux across any remote surface, Sqo, surrounding 
a volume Voo containing and moving with the solitary wave, it is, to second order in ij), the 
same as the mass flux implied by ipQ. This implies that 

I [$*V^-4>V^*) -2u7|^| 2 £]-dS = 0. (22) 

This in turn implies that the operator L is hermitian: 



ipL*i>*dV = / i)*Li)dV. (23) 

From this relation and ()15|) . it follows that 

/ [\Lp\ 2 - ^lp*L*p*}dV = -2ia [ P L*p*dV 

= -2i(T f p*LpdV=-2ia [ p*[^lp* -2\ap]dV. (24) 

Re-arranging this equation and using (|15|) again, we find that 

/ [\Lp\ 2 -|Vo| 4 M 2 ]^ = -Aa 2 j \p\ 2 dV. (25) 

It follows that a 2 can only take real values. The numerical results indicate that the eigen- 
value spectrum is infinite and discrete and (|25jl shows that its limit point is at a 2 = — oo. 
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But only the necessarily finite number of positive a 2 are relevant to the physical problem; 
see §2.3. 

For the upper branch of the JR dispersion curve, we found only one family of con- 
vincingly positive a 2 , but this family (for which m = 0) is sufficient to establish that the 
rarefaction waves are unstable to axisymmetric perturbations. The family seemed to dis- 
appear as the cusp U = U c was approached, and to become continuous with a family of 
negative a 2 modes on the lower branch of the JR dispersion curve. Numerical work is inca- 
pable of establishing this fact unequivocally but an analytic argument given in Appendix A 
lends support. As U — > c, the growth rate of the instability tends to zero too (see Appendix 
B). One of the more striking results of the numerical work was the discovery that a is small 
everywhere on the upper branch. It was found that a has a single maximum a max (U), of 
approximately 0.012, attained for U ~ 0.68. The fact that a max (U) is so small means that 
the enfolding time on which the instability grows is long, being greater than 1/0.012 ~ 84 
in all cases. Fig. |3] shows the density and the phase contour plots of the wave function for 
the the fastest growing mode of perturbation for rarefaction solitary wave moving in the 
positive z— direction with the velocity U = 0.68. Fig. 2] depicts the maximum growth rate 
<T max ([7) as a function of the rarefaction wave velocity, U. 

No instabilities were discovered for the asymmetric (m 7^ 0) modes on either branch 
of the JR dispersion curve, but the neutral \m\ = 1 modes, corresponding to infinitesimal 
displacement of the basic state in a direction perpendicular to Oz, were located. 

2.3 Complex growth rates 

When a is complex, the complex conjugate of either of (|11|) contradicts the other. On 
making the exchange 1 <-> 2 and taking the complex conjugates of both equations (|12j) . we 
see however that, if a is an eigenvalue, so is a* . Thus solutions to @ exist of the form 

£(x,t) =pi(x)e CT *+p 2 (x)e CT **, ^*(x,t) = ^(x)e CT * + p 1 (x)e <T ** . (26) 

By making use of the Hermitian property of L, it follows from an argument similar to 
the one that led to (|25|) that 

J\\L Pl \ 2 -^\ Pl \ 2 -i; 2 (p* 2 L*pl-plL*p* 2 )]dV = -4a 2 j \ Pl \ 2 dV , (27) 

J [\Lp 2 \ 2 - |V>o| 4 |p 2 | 2 - r 2 (PiL P 2 -p 2 L Pl )]dV = -4a 2 J \p 2 \ 2 dV . (28) 

On subtracting corresponding sides of these equations we obtain 

4^ 2 A|p 2 | 2 - IpiI 2 ]^ = J [\L Pl \ 2 - \l P2 \ 2 - |Vol 4 (M 2 - N 2 ) 

+ r 2 (PiL P 2 ~ P2L Pl ) + ^l(jplL*pl - P * 2 L*pl)]dV . (29) 



8 



Figure 3: The density (a) and the phase (b) of the wave function for the the fastest growing 
mode of perturbation for rarefaction solitary wave moving in the positive z— direction with 
the velocity U = 0.68. (These results were derived by solving the spectral problem (|16|) - 
(|17|) numerically). 
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In what follows, we shall assume that there are no eigenvalues (apart from the real a 
of §2.1) for which 



/ 



\ Pl \ 2 dV = J \ P2 \ 2 dV, (30) 

[|L m | 2 -|Lp 2 | 2 -^ | 4 (bl| 2 - \P2\ 2 ) 

+% 2 ( Vl Lp 2 -p 2 L Pl ) + i>l{ V \L*p* 2 - P * 2 L* V \)]dV = Q. (31) 



hold simultaneously. It then follows from ([29j) that a 2 is real. Therefore the only eigenvalues 
that are complex are purely imaginary. The main import of this result is that §2.2 has 
already located all modes that can become unstable. 
For modes with imaginary a = iu, we substitute 

Pl = ui + iv 1 , p 2 = u 2 + iv 2 , (32) 

into ([12[) to obtain an 8 th order system that replaces the 4 th order system ([16[) and ([170: 

(33) 
(34) 
(35) 
(36) 



where 



2.4 Large r behaviour 
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+ 2Ud z vi - 


iul- 


vl)u 2 
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— 2uoVou 2 = 
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{ul- 
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—2oju 2 
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- 2Ud z u 2 + 
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— 2uqVqU\ = 


—2ujv 2 . 
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= V 2 


+ 1- 


2u\ - 2v\ . 





(37) 



It might seem that (jlOjl is the obvious and correct boundary condition to apply to solutions 
of ©, but there are some subtleties. Consider first the solutions developed in §2.3. For 
r — > oo, «i and V\ make contributions to ip that are dominantly proportional to 

r -1 exp[i(/cr + at)] ■ (38) 

Since uq — 1 = 0(r -3 ) and v$ = 0(r~ 2 ), the possible values of k implied by ([3*51 - l[36[) 
are: 

oj±wk = ±\k^J{k 2 + 2), (39) 

where w = Uz/r (—l/y/2 < w < l/y/2). The eigenfunction is a linear combination 
therefore of 8 solutions of the form ((38|) , for the 8 roots of the 2 quartic equations 

k 4 + 2(1 - 2w 2 )k 2 ± 8u>wk - 4w 2 = . (40) 
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It is easily shown that 4 of these roots are real; 2 are positive and 2 negative. These 
correspond to ingoing and outgoing sound waves. They satisfy the requirement (|1U|). but 
they do not obey the condition (|22[) needed to establish that L is self-adjoint. This difficulty 
can be surmounted however by requiring that the waves are reflected from Sqq. A more 
serious issue concerns the remaining 4 complex roots of (|4(Jj) . two of which have positive 
real parts and two negative. To satisfy (|lUj) . the eigenfunction must exclude the two roots 
for which Im(/c) < 0, irrespective of whether the remaining k correspond to ingoing or 
outgoing energy flux. 

In short, hidden in the succinct statement (|l()j) . are physical requirements on the eigen- 
functions that are far from obvious and may be controversial. The same is true in the 
case of real a considered in §2.2. For this reason, and also to throw further light on the 
evolution of perturbations to the basic state, it was thought desirable to carry out the nu- 
merical experiments reported in §3 below. These attack the initial value problem without 
preconceptions such as (10), and without the necessity of supposing that the perturbation 
to ipo is infinitesimal. 

3 Numerical simulations involving rarefaction waves 

In this section we shall report on numerical solutions of the GP equation (JIJ that show 
that the solitary waves on the upper branch can, despite their instability, be surprisingly 
robust, but that nevertheless they eventually evolve onto the lower branch. 

3.1 Perturbed unstable solitary wave 

Starting with a rarefaction solitary wave from the upper branch of the JP dispersion curve, 
we consider perturbations of the form 

ip(x, t = 0) = Vo(x) + AA(x), where N = ^ a k exp(zk • x) . (41) 

k 

Here J\f corresponds to the state of a weakly interacting Bose gas in the kinetic regime 
jH]; the phases of the complex amplitudes ak are distributed randomly. To make the 
disturbance decay to zero far from the solitary wave and to give it both axisymmetric and 
non-axisymmetric components, we take 

1 M 

k=(i,m,n) 

x exp[-0.01(x 2 + 2y 2 + 3z 2 )]. (42) 

The maximum frequency, M, and the modulus of the complex amplitude |ok| are two 
parameters that we vary. 
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In what follows we say that a perturbed solution evolves to a JR solitary wave of 
velocity U if, in a sufficiently large neighbourhood of its centre and for a significant period 
of time, it becomes and remains axisymmetric and form-preserving to good accuracy, and 
if in addition its minimum u corresponds to the U shown in Fig. [3 By "a significant 
period", we mean an elapsed time exceeding 10, by "good accuracy", we mean that the 
density deviations do not exceed 0.01 and, by "a sufficiently large neighbourhood", we 
mean within a sphere of radius 20 surrounding the origin. 

We performed numerical integrations of ()20j) with the time derivative —2\dt^ restored 
to the left-hand side. Since the integration was terminated at a finite time t, before sound 
waves from the disturbance can reach r = oo, the boundary condition (|1U|) on ip is irrelevant. 
We employed the same mapping as in §2.2. We used fourth order finite differences in space 
except for second order finite differences at the boundary, and fourth order Runge-Kutta 
integration in time. The initial state at t = was given by (|41|) - (|42[) . where V'o( x ) are 
rarefaction waves from the upper branch of the JR dispersion curve. We principally focused 
on two cases: U = 0.68 and U = 0.69; see Fig. ^ 

In each case, our simulations show that, for 5 < M < 20 and \a^\ < 4, the wave has, by 
t ~ 20, rid itself of the perturbation, apparently by acoustic radiation, and has returned 
to the vicinity of the upper branch, though with a velocity that is slightly smaller (by 
1% or less) than that of the initial tpQ in (|41j) . Fig. [5] shows snapshots of the evolution 
of the rarefaction waves for U = 0.68 and U = 0.69, for the initial perturbation (|41|) 
- (|42j) with M = 20 and |ak| = 2. By t = 18, each solution had evolved close to an 
unperturbed rarefaction wave, and at t = 36 it still remains in its vicinity. On a longer 
time scale (t > 100), the rarefaction waves either lose their energy and momentum to sound 
or collapse onto the lower branch of the JR dispersion curve. 



3.2 Evolution of unstable rarefaction waves into vortex rings 

In [7] an algorithm was developed for generalised rational approximations to the JR solitary 
waves having the correct asymptotic behaviour at infinity. An axisymmetric solitary wave 
moving with velocity U along the z— axis is well approximated by ip(s, z) = 1 + u(s, z) + 
iv(s,z) where 

aoo + aioz' 2 + a 0l s 2 + mc 2 7 / 4 U(2z' 2 - f(U)s 2 )(z' 2 + f{U)s 2 ) 
(1 + c w z' 2 + c 01 s 2 + c 20 (z' 2 + /(C/> 2 ) 2 ) 7 /4 

, &oo + b w z' 2 + b 01 s 2 - mc^V 2 + f(U)s 2 ) 2 
2 (l + c 10 z' 2 + c 01 s 2 + c 20 (z' 2 + f(U)s 2 ) 2 y^ [ } 

where z' = z — Ut and f(U) = 1 — 2U 2 . Here dy, bij, Cij and the dipole moment m 
are functions of U determined by series expansions. For example, the rarefaction wave 
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Figure 5: (Colour online) Snapshots at t = (black line), t = 18 (dark grey or blue line) 
and t = 36 (light grey or red line) of the density of the condensate along the direction of 
the motion. Here |^>(0, 0, z)\ 2 was obtained by numerically integrating the GP model (pQ) 
for the perturbed rarefaction waves. 
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U = 0.69 is given by 



-0.2779-0.00182z /2 -0.00128s 2 +0.00035(2z /2 -0.0478s 2 )(z /2 +0.0478s 2 



-0.34761 -0.02198z 2 - 0.00262s 2 -0.00051 



V 2 + 0.0478s 2 ) 2 



V=(l + 0.11749 z' 2 + 0.01470 s 2 + 0.00356(z' 2 + 0.0478 s 2 ) 2 ) 



where 

By © and the momentum and energy are 
p = 2-ir / (uv x — vu x )s ds dx , 



7/4 



r 2 2 , 2 , 2 



i(2u + u 2 + t> 2 ) 2 ]sds<ix, 



(44) 
(45) 

(46) 
(47) 



where u z = du/dz, u s = du/ds, etc. By substituting the approximation (|4*4*|) into (fl6|) and 
(|47|) . we obtain p ~ 84.15 and £ ~ 61.2; these may be compared with the values p ~ 83.2 
and £ ~ 60.0 obtained by from direct integration of ©. The maximum residual error, 
calculated as the global maximum of the square of the error amplitude [see (29) of [7]] is 
~ 0.003. Evidently, (fi4*|) provides a starting point that is very close to the upper branch. 

To follow the evolution of (J44[) over long times, we performed a numerical integration of 
(^Q), using a finite differences scheme with open boundaries to allow sound waves to escape 
(see our previous work (Sj or (Hj for a detailed description of the numerical method). The 
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computational box has dimensions D = 100 , the space discretization was dx = dy = 
dz = 0.5, and the time step was dt = 0.1/ (dx 2 + dy 2 + dz 2 ). The small initial deviation of 
(|44|) from the upper branch grows with time, and the solution slowly evolves into a vortex 
ring. Fig. illustrates this process. Circulation is acquired at t ~ 148 and a well-formed 
vortex ring emerges by t ~ 160. 

Figure 6: (Colour online) Snapshots of the contour plots of the density cross-section of a 
condensate, obtained by numerically integrating the GP model (^Q). The initial condition 
was ip = 1 + u + iv, with u and v given by 1)44(1 . Black dashed lines show zeros of the 
real and imaginary parts of tp at t > 0. Their intersection therefore shows the position of 
topological zeros. Both low and high density regions are shown in darker shades so that 
regions of intermediate density are emphasised. Only a portion of an actual computational 
box is shown. 




t t 90 




t = 148 t = 166 



3.3 Vortex nucleation resulting from the interaction of rarefaction waves 

In |7] we considered the evolution of two identical rarefaction waves (U = 0.63) from the 
lower branch. These were centred on the z— axis, and moved along it, one behind the 
other. It was found that the trailing wave transfered part of its energy and momentum 
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to the leading wave, until eventually the latter became a vortex ring. We now consider 
an analogous situation but one in which both rarefaction waves initially lie on the upper 
branch. We numerically integrated Q in the frame of reference moving with U = 0.69, 
taking the initial state to be ip(t = 0) = tpo(x — 5, s)tpo(x + 5, s), where tpo(x, s) is defined 
by (|44p. Fig. [7| shows snapshots of the subsequent evolution. This is very different from 
the evolution of the two rarefaction waves from the lower branch, although there is again a 
transfer of energy and momentum from the trailing wave to the leading wave. The leading 
wave does not, however, evolve to a higher energy state on the upper branch but instead 
slowly transforms itself into a vortex ring on the lower branch. Apparently, the transfer of 
energy and momentum from the trailing vortex is insufficient to carry the wave to a higher 
energy state on the upper branch; it remains 'beneath' that branch and can therefore only 
evolve onto the lower branch; see Fig. 

3.4 Evolution of the density depletions in condensates 

There are several ways of creating rarefaction waves in a condensate, that make use of some 
of the methods described in |lUj . A tangle of vortices is created when angular momentum is 
transmitted to the condensate by rotationally stirring it with a laser beam [IT . In addition 
to vortices, such stirring creates many local depletions in the condensate density. It was 
shown in ^U] that such depletions are unstable and generate both vortices and rarefaction 
waves. We shall now show that, by creating a shallow ellipsoidal depletion in condensate 
density, it is possible to generate rarefaction waves alone on the upper branch of the JR 
dispersion curve . 

We numerically integrated (^Q) taking as our starting point a spheroidal depletion of 
condensate: 

V>(x, 0) = \ + \ tanh [0.01(z 2 + 0.5s 2 - 36)] . (48) 

The minimum amplitude of the initial state is 0.33, at the origin. After the condensate 
fills the cavity, it begins to expand, its density oscillating around the unperturbed state 
ifj = 1. These depressions in density are unstable, in a manner similar to the instabil- 
ity of Kadomtsev-Petviashvili 2D solitons in 3D [12|. This results in the creation of two 
rarefaction waves moving outwards in opposite directions. Fig. El shows snapshots of the 
density in the z = cross-section. Well-defined vorticity-free localised axisymmetric dis- 
turbances - rarefaction waves — are formed at around t ~ 30. Their axes of symmetry lie 
on Im(ip) = v (x, t) = and are plotted as solid lines in the t = 33 snapshot. 

Various stages in the collapse of a stationary spherically symmetric bubble were elu- 
cidated in |T2], where conditions necessary for vortex nucleation were also established. A 
similar analysis applies to the evolution of the density depletion (|48ft . for which the density 
is everywhere nonzero initially. In particular, it would be possible to calculate the critical 
radius of a spherical depletion, as a function of minimum density, for rarefaction waves to 
be nucleated. 
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Figure 7: (Colour online) Snapshots of the contour plots of the density cross-section of a 
condensate, obtained by numerically integrating the GP model (^Q). The initial condition 
was ip(t = 0) = ipo(x — 5, s)ipo (x + 5, s), given by (fH)) where tpo = 1 + u + iv. Black dashed 
lines show zeros of the real and imaginary parts of ip. Their intersection therefore shows 
the position of topological zeros. Both low and high density regions are shown in darker 
shades so that regions of intermediate density are emphasised. Only a portion of an actual 
computational box is shown. 




t = 75 t = 90 



4 Conclusions 

We have analysed the linear stability problem for solitary waves of the GP equation. We 
have shown that the rarefaction solitary waves on the upper branch of the JR dispersion 
curve are unstable to infinitesimal perturbations and we have calculated the maximum 
growth rate, and found it to be small, b eing approximately <r max — V2 x 0.012c/£ ps 70 
s _1 (taking the healing length £ = 0.7//m and the Bogoliubov speed of sound c = 2.8 
mm s _1 as in NIST experiments 2J. Our results indicate that the solitary waves on the 
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Figure 8: (Colour online) Snapshots of the contour plots of the density cross-section of 
a condensate obtained by numerically integrating the GP model (0) starting from the 
spheroidal density depletion (|48|) . Black solid lines show zeros imaginary parts of ip, the 
straight segment being the axis of symmetry of the rarefaction wave. Both low and high 
density regions are shown in darker shades so that regions of intermediate density are 
emphasised. 




f V> 




lower branch and all two-dimensional solitary waves are linearly stable. By numerically 
integrating the GP equation, we studied the properties of the unstable solitary waves and, 
in particular, we showed that the system tends to remain close to the unstable solution, 
spending a significant amount of time, of the order of a few multiples of r m ; n = l/<r max , in 
its vicinity. Ultimately however it collapsed onto the lower branch or broke up into sound 
waves. 

These conclusions suggest that rarefaction solitary waves could be experimentally de- 
tected and and studied in Bose condensates alongside the vortex rings. They can be created 
by laser beams and evolve from a local depletion of the condensate. Rarefaction waves will 
also appear as the result of the self-evolution of a strongly nonequilibrated Bose gas |14j . 
Like the vortex rings, they can in principle be observed in expanding condensates. 
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Appendix A: Eigenvalues at the cusp 



It was argued in §2.1 that, for every U, the spectrum of a 2 , where a is the growth rate of 
perturbations is real, discrete and infinite, with limit point at a 2 = — oo, i.e., only a finite 
number of positive a 2 can exist for each m. The largest of these is the most significant 
since, when positive, it gives the largest growth rate for instability. Our numerical work 
strongly suggested that this a 2 is part of the spectrum for the axisymmetric mode and that 
it changes sign precisely at the cusp, U = U c , being positive for U > U c and negative for 
U < U c . The aim of this Appendix is to re-enforce this conclusion by showing analytically 
that this a 2 can change sign only at the cusp. The argument is analogous to one presented 
by Kuznetsov and Rasmussen ^7] for 2D solitary waves. 

In the axisymmetric case (m = 0), (|15j) always possesses a neutral solution (a = 0) 
corresponding to an infinitesimal displacement of tpo in the z— direction; we write this as 

= d?p /dz . (49) 

We now seek a neighbouring solution for which a is small, vanishing when a double root 
of a exists. We expand this solution as 

p = p W +(Tp W +a 2 p (2) + _^ . (50) 

Substituting (J5U|) into (|T5|) . we obtain 

Lp(°)-V>oV 0) * = 0, (51) 
Lp W - ?PoP {1> = -2ip(°\ (52) 
Lp( 2 )-VoV 2) * = -2fp {1) . (53) 
The differential of ()15|) with respect to U and (|49|) show that the solution to (|52j) is 

p« = -dipo/dU. (54) 

A consistency condition follows from (J53|). It is easily seen, from ()13|) and (J51|) and the 
condition (jlUj) on tp at infinity, that 

p(0)*(£ p (2) _ ^(2)*) + p (0) {L * p (2h _ ^(2))] d y = . (55) 



It now follows from (J53*|) that 

dp _ 1 f ( dipQ dif) dip dipt 



dU i / V dz dU dz dU 



dV = 0. (56) 



The left-hand equality here is established by differentiating © with respect to U and 
carrying out integrations by parts. 

It follows from (|56|) that the double zero of a can occur only at the cusp, the only point 
on the JR dispersion curve where p'(U) is zero. 
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Appendix B: The Gross-Pitaevskii equation for U — > c 

In multidimensions the density depletion of the solitary rarefaction wave as well as the 
characteristic inverse length along the axis of propagation tend to zero as the solitary wave 
velocity, U, approaches the speed of sound, c. This allows us to reduce the GP equation 
to the Kadomtsev-Petviashvili type I (KPI) equation, that describes the propagation of 
acoustic waves of small amplitude and positive dispersion. 

In [3] such a reduction was obtained as a compatibility condition of the equations 
written to 0{e ) where the small parameter e ~ y/(c — U). In what follows we derive the 
eigenvalue problem of the GP equation in the limit of U — > c and show that it coincides with 
the eigenvalue problem for the KPI equation for an appropriate scaling of the eigenvalues. 

We start by separating the real and imaginary parts of the basic state as V>o = u o + i^o 
and of the perturbation as ip = u + \v; see Section 2. The basic solution satisfies © so 
that 

V 2 u + (1 - ul - v$)u + 2Ud z v = 0, V 2 v + (1 - tig - ug)«o - 2Ud z u = 0, (57) 

while u and v satisfy (|16|) - (|17j) where, without loss of generality, we take m = 0. To make 
the reduction to the KPI equation we write 



z/e, 



so that 



d z 



ed z , V s 



e 2 a 2 z + e 4 v 2 



Hi 



(58) 



where V 2 ^ 



d 2 + s 1 d s . Next we make the transformation 



u 



1 + e u , 



eu, 



U^c + e 2 U, a 



eV. 



(59) 



Note that the transformation of the eigenvalue a defines a slow timescale, the only one that 
makes the time-dependent GP equation consistent with the time-dependent KPI equation; 
see below. 

The equations governing the basic solution and the perturbations become 



1 e 2 
d z v - V2u - -j=vl + -= 



d 2 u Q - (3u + vfyuo + 2Ud z v 



d 2 v - V2d z u - (2u + vl)v + e 2 



V 2 H v - ulvo - 2Ud z u 



0, 

: 0, 



2{u + vov) + —= 
\J2 



d z u — (6uq + Vq)u — 2uqvqv + 2Ud z v — 2av 



(60) 
(61) 



+ 0(e 4 ) K62) 



d 2 v — V2d z u — (2uq + 3vq)v — 2vqu 



V\v — UqV — 2uoVqu — 2U d z u + 2au 



+ 0(e 4 ) = 0. 



(63) 
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Now we expand all functions and the velocity U in powers of e 2 as uq = Uq + e 2 ^o, etc. 
At leading order, (|S0|l - give 



d z v? 



■ z u - V2u° - — = 0, 



v/2 



d 2 v° - V2d z u° - (2u° + v° 2 )v° = 0, 

d z v° - V2u° - V2v%v° = 0, 

d 2 z v° - V2d z u° - (2u° + 3v° 2 )v° - 2v° u° = 0. 

We first consider (|6~H) - (ffil))) and notice that implies (JHSJ), so that 



i/° - — -r) 7>° 



,,02 



V2 



(64) 

(65) 
(66) 
(67) 



(68) 



We use (j68fl to eliminate Uq from the leading order expressions in square brackets in 1)60(1 
and 1)61)1 which we denote by R and S respectively. After making simplifications we obtain 



R = d 2 u° - (3u° + v$ 2 )u° + 2U°d z v° 



°z v v 0°z v 



S 



1 

V | Uq - u° 2 v ° - 2U°d z u° 



2^0 ' 



0q „,0 



Equations (JHUj) - (|61)1 become 



d z vl - V2ul - V2v%vl = — r R, 

V2 

dlv l - V2d z ul - {2u° + 3^)^ - 2v Q Q ul = -S. 
The compatibility of (|71|) and ()72j) implies 

5 = -^d 2 i? + v%R, 
v2 

which is the KPI equation as derived in 



0. 



(69) 
(70) 

(71) 
(72) 

(73) 
(74) 



Next we consider (|62|) and (|63|) . At the leading order both (|66|) and ()67)) require that 

1 

72 



- 9 ..0..0 



Voir 



(75) 
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Note that, if we had assumed a different scaling for a, we would now face an inconsistency. 
For example, if we had adopted a = 0(e), we would have found at this point in the 
argument that a = 0, i.e., that a is asymptotically smaller than the assumed 0(e). 

At the next order beyond and (|67|). expressions R' and S' analogous to R and S 
arise: 

R> = ^ _ ( 6n + v 02 )u _ 2u v v + 2[/ a ^0 _ 2(JV + (?6) 

5' = V^u - u§V - 2u° v%u° - 2U°d z u° + 2au° + (2u\ + 6^)v° + 2^n°. (77) 

We use (fo%)) and ((751) to eliminate u[] and u° from (|76() and (|77j). The equations and 
(|67l) then have the form 

d^ 1 - + wgv 1 ) = -^=R', (78) 

v2 

flgt; 1 - V28.U 1 - (2u° + 3^ 2 y - 2v%u 1 = -S', (79) 
and the compatibility of (|78[) and (f?9"|) leads to the eigenvalue problem 

S' = -^d z R' + v° R', (80) 
v2 

or, after simplification, 

2v / 2t/ dV - sV^dzidzV^dzV 1 ] + -figy 1 - V^u 1 = 2\[2od z v x . (81) 

Notice, the spectral problem (|81|) is the same one as obtained from (|74j) written in time- 
dependent form as 



Z 



2V23 t V - 2V2U°d z V + -^(d z V) 2 - \dlV 

v2 ^ 



+ V^y = 0, (82) 



where V = v^ + v 1 exp(at). This proves that the KPI equation and the GP equation share 
the same linear stability properties, the eigenvalues being related by ctkp(c — Ugp) 3 ^ 2 = 
a GP as U GP -> c. 
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